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LONG-TERM  GOALS 

1.  Application  of  the  inverse  scattering  transform  (1ST)  to  the  time  series  analysis  of  laboratory  and 
oceanic  wave  data.  The  approach  may  be  viewed  as  a  generalization  of  linear  Fourier  analysis  and 
is  loosely  referred  to  as  "Nonlinear  Fourier  Analysis  or  Generalized  Fourier  Transform"  (GFT). 

2.  A  major  focus  has  been  the  application  of  1ST  to  the  study  of  “rogue,  freak  or  giant”  ocean  waves. 
The  emphasis  has  been  on  the  study  of  the  physical  mechanisms  leading  to  the  generation  of  rogue 
waves  in  random  sea  states. 

3.  A  third  long  term  goal  is  the  development  of fast  algorithms  for  numerically  integrating  the 
space/time  dynamics  of  deep-water  wave  trains.  While  1ST  is  limited  to  the  numerical/analytic 
integration  of  the  so-called  “soliton  equations,”  I  have  discovered  how  the  GFT  can  be  used  to  solved 
higher  order  equations  for  which  study  of  the  dynamics  have  previously  been  limited  to  numerical  (as 
opposed  to  analytical)  approaches.  I  discuss  herein  how  the  GFT  can  be  used  both  for  the  analytical 
study  and  extremely  fast  numerical  integration  of  the  extended  nonlinear  Schroedinger  equation  for 
fully  three  dimensional  wave  motion.  I  discuss  a  number  of  considerations  for  the  development  of 
extremely  fast  GFT  algorithms,  which  I  refer  to  as  the  xGFT. 

OBJECTIVES 

1.  The  objective  of  the  present  research  program  is  the  development  of  fast  numerical 
multidimensional  Fourier  techniques  applied  to  a  wide  range  of  wave  modeling  and  analysis 
problems. 

2.  Important  progress  made  in  the  past  year  has  been  the  development  of  new  algorithms  for 
multidimensional  Fourier  analysis.  These  algorithms  are  the  key  to  future  hyperfast  applications  of 
the  method. 

APPROACH 

The  nonlinear  modeling  project  that  I  discuss  herein  is  quite  new  and,  remarkably,  has  an  overall 
structure  as  seen  in  Fig.  1.  Basically  this  means  that  the  first  step  consists  of  a  preprocessor  part  which 
determines  the  nonlinear  time  evolution  of  the  linear  Fourier  spectrum,  i.e.  it  computes  the  linear 
Fourier  spectrum  at  some  large  number  (say  N  ~  1000 )  0f  desired  time  steps  (one  normally  takes  a 
sufficient  number  of  time  steps  to  make  a  movie  for  studying  the  nonlinear  wave  behavior).  The 
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second  step  is  to  take  the  inverse  FFT  of  the  time  varying  Fourier  transforms.  Additionally  the 
mathematics  (soliton  theory,  Riemann  theta  functions  and  algebraic  geometry)  is  not  that  ordinarily 
studied  by  physical  oceanographers  and  this  may  be  an  impediment  that  could  test  the  patience  of 
those  who  might  want  to  apply  the  method  (see  for  example  the  book  by  [Baker,  1897]).  As  we  shall 
see  however,  the  speed  gains  are  quite  dramatic,  i.e.  two  or  three  orders  of  magnitude  faster  than  the 
direct  numerical  integration  of  nonlinear  wave  equations  by  the  FFT.  The  preprocessor  is  unique  for 
each  integrable  wave  equation  as  discussed  in  detail  below.  Furthermore,  for  each  nonintegrable 
equation  a  sub-processor  must  be  added  as  shown  in  the  blue  box  of  Fig.  1. 

We  first  consider  the  Kadomtsev-Petviashvili  (KP)  equation 

Vt  +  corlx  +  «»77x  +  PVxxx  +  Y  ^xdyy  =  0  ^ 

Here  7CATT)  is  the  wave  amplitude  as  a  function  of  the  two  spatial  variables,  *,  -v  and  time,  t.  The 
KP  equation  (1)  is  a  natural  two-space-dimension  extension  of  the  KdV  equation.  The  periodic  KP 
solutions  include  directional  spreading  in  the  wave  field: 


q(x,t)  =  2 


d2 

— j\nO(x,y,t  |B,K,Q,0) 
dx 1 


(2) 


Here  the  generalized  Fourier  series  has  the  form  given  in  (4)  below,  where  the  phase  has  the  two 
dimensional  expression: 


X(x,y,t)  =  Kx  +  Ly-  Clt  -  O 


(3) 


Here  the  spatial  term  K*  has  been  joined  by  the  lateral  spatial  term  which  allows  wave  spreading 
to  be  taken  into  account.  The  KP  equation  is  the  first  nonlinear  step  toward  a  directional  sea  state;  KP  is 
however  limited  to  small  directional  spreading.  Improving  the  directional  spreading  characteristics  of 
the  KP  equation  requires  the  addition  of  physically  important  corrections  to  the  equation:  I  give 
numerical  examples  below. 

The  generalized  Fourier  series,  Pxd  |B,k,<o,<|>)  ^  js  gjven  by  the  expression 


aaj  y^j 

0(x,y,t  I  B,<j>)=  ^  X  -  X 


N  Y  N  N 

i^mnXn X  X  mm  Mn 

?  «= i  2  m=\  „=i 


m,  — — go  — — go  m  =-oo 


(4) 


where  ^"x  +  0)id  +  fin  .  The  function  ^xd’d  |B,(j))  js  a|so  capec|  a  Riemann  theta  function 
or  multidimensional  Fourier  series.  Here  B  is  the  Riemann  matrix  (the  “spectrum”  of  the  solution),  the 

vectors  k,  1  constitutes  the  usual  wave  numbers,  the  vector  ®  contains  the  frequencies  and  the  vector  ^ 
forms  the  phases.  The  inverse  problem  associated  with  (2),  (3)  allows  one  to  determine  the  Riemann 
matrix,  wave  numbers,  frequencies  and  phases  appropriate  for  solving  the  Cauchy  problem  for  KdV: 

Given  the  spatial  variation  of  the  solutions  at  { =  0,  ?X-Ty,0)  ^  compUte  the  solution  for  all  time, 

>l(x->yd)  This  is  a  necessary  step  for  the  numerical  simulations  presented  heein.  The  solitons,  Stokes 
waves  and  sine  waves  lie  on  the  diagonal  of  the  Riemann  matrix ;  the  off-diagonal  terms  contain  the 
nonlinear  interactions. 
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Why  is  the  above  approach  useful  for  hyperfast  numerical  simulations?  Because  the  Riemann  theta 
function  can  be  programmed  as  a  fast  theta  function  transform  (FTFT),  just  as  the  Fourier  transform 
can  be  programmed  as  a  fast  Fourier  transform  (FFT).  Therefore  the  numerical  integration  of  KdV  (1) 
can  be  evaluated  at  specific  time  points,  necessary  only  for  graphical  purposes  or  for  extracting  useful 
properties  of  the  sea  surface.  This  contrasts  to  the  FFT  that  must  be  evaluated  at  very  small  time  steps 
when  used  for  the  numerical  integration  of  a  nonlinear  partial  differential  equation. 

The  unidirectional,  deep-water  case  is  governed  by  the  nonlinear  Schroedinger  equation : 

Kv,+Cgyx)+fiy/„  +  v  HV=0  ro 


where  the  coefficients  are  computed  in  the  usual  way  in  terms  of  the  carrier  wave  number  and 
frequency.  The  spectral  solution  to  this  equation  is  given  by 


yKx,t) 


^0(x,t  | 

0(x,t  |  B,8+) 


(6) 


Thus  the  solution  is  the  ratio  of  two  multi-dimensional  Fourier  series  with  different  phases, 
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Solving  Nonlinear  PDEs  With 
The  Multi-dimensional  Fourier 
Transform 


Fig.  1  Simple  flow  chart  of  the  proposed  numerical  models  for  nonlinear 
partial  differential  equations  with  a  well-defined  dispersion  relation. 
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For  directionally  spread  waves  in  deep-water  the  nonlinear  Schroedinger  equation  has  the  form: 


K ¥,+Cgif/x)  +  juif/xx  -  2jU<//yv  +  v\  iy\ 1  y/=  0 


(7) 


An  important  result,  fundamental  in  the  present  research  program,  is  that  the  above  relatively  simple 
cases  form  building  blocks  for  the  solution  to  more  complex  equations  such  as  the  Boussinesq 
equations  (see  Fig.  1)  and  the  full  Euler  equations  (see  below).  Another  important  aspect,  to  be 
addressed  in  the  coming  year,  is  that  variable  shaped  boundaries,  variable  bathymetry  and  wind  and 
dissipation  forces  can  be  easily  added  to  the  above  equations  and  to  the  Euler  equations. 

WORK  COMPLETED 

Up  to  this  point  I  have  not  given  any  of  the  details  of  the  nonlinear  preprocessor  in  the  new  approach 
for  obtaining  numerical  solutions  to  nonlinear  wave  equations.  I  will  now  do  so,  but  without  many  of 
the  mathematical  details.  First  note  that  one  begins  with  the  Riemann  theta  functions  and  makes  a 
transformation.  One  of  the  simplest  types  of  transformation  is  given  by  (2).  Second  I  will  use  the 
following  result,  i.e.  that  the  theta  function  can  be  reduced  to  an  ordinary  linear  Fourier  series  with  time 
varying  coefficients: 

OO  00 

(Kx,y,t)=  £  Z 

tyi= — oo  n= — GO 


Here  the  time  varying  coefficients  can  be  written  analytically  in  terms  of  the  Riemann  matrix 

and  phases.  Note  that  the  above  expression  is  rather  simple,  i.e.  it  is  just  an  ordinary  linear  Fourier 

series,  which  in  numerical  applications  has  ~  N  terms.  The  theta  function  itself  is  quite  different  in 

behavior  for  it  has  an  exponential  number  of  terms  ~  10^ ,  which  can  be  a  huge  number.  Consider  for 
example  a  case  where  N  =  30  s  i.e.  30  cnoidal  waves  in  the  spectrum.  Then  the  number  of  terms  in  the 
theta  function  is  far  greater  then  the  number  of  grains  of  sand  on  the  earth,  the  number  of  galaxies  in 
the  universe,  greater  than  Avogardro’s  number  and  greater  then  the  number  of  seconds  since  the  big 

bang!  Clearly  the  preprocessor  has  got  to  be  pretty  efficient  to  reduce  10^  (~1030 )  to  ~N2  (900)!  I 
call  the  mathematical  process  to  reduce  the  theta  function  onto  the  ordinary  linear  Fourier  modes  a 
“collapse”  in  analogy  with  quantum  physics  when  one  refers  to  the  collapse  of  a  wave  function  during 
the  process  of  a  physical  measurement.  Indeed  we  say  that  we  collapse  the  theta  function  onto  the 
linear  Fourier  modes.  This  process  requires  some  knowledge  of  the  algebraic  geometry  of  theta 
functions  [Baker,  1 897]  that  I  do  not  go  into  here. 

I  now  consider  some  interesting  numerical  results.  To  be  concrete  I  have  chosen  the  KP  equation,  i.e. 
directionally  spread  waves  in  shallow  water  that  are  distributed  over  relatively  small  angles.  An 
example  of  a  directionally  spread  wave  train  is  given  in  Fig.  2.  The  significant  wave  height  is  1.5  m  and 
the  water  depth  is  8  m.  There  are  20  cnoidal  waves  in  the  spectrum  with  a  maximum  modulus  of  0.84 
(these  are  strongly  nonlinear  waves,  but  they  are  still  less  nonlinear  than  solitons  which  have  a  modulus 
of  1).  Because  the  moduli  of  the  spectral  components  is  relatively  large  the  waves  are  really  Stokes-like 
in  their  shape,  a  result  of  the  nonlinearities  in  the  KP  equation. 
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I  now  turn  to  the  modified  or  extended  KP  equation,  which  I  refer  to  as  xKP.  This  equation  is  found 
from  the  KP  equation  by  extending  it  an  additional  two  orders  of  approximation  beyond  the  KP 
equation  and  for  many  purposes  xKP  behaves  essentially  like  the  Boussinesq  equation.  To  add  this 
extended  feature  to  the  model  in  Fig.  1  we  must  of  course  add  a  time  varying  Riemann  matrix  and 
phases.  To  this  end  the  model  of  Fig.  1  must  be  modified  by  adding  the  dynamics  of  the  Riemann 
matrix  and  phases  to  the  blue  box. 

The  numerical  results  for  the  Boussinesq  model  are  shown  in  Figs.  3  and  4.  One  of  the  major  new 
results  which  have  come  from  these  simulations  is  the  appearance  of  a  new  kind  of  rogue  wave  which 
occurs  only  in  shallow  water  wave  dynamics  and  is  not  related  to  the  Benjamin-Feir  instability. 

Shallow  water  rogue  waves  of  the  type  observed  in  the  simulations  occur  due  to  strong  nonlinear 
interactions  between  two  or  more  cnoidal  waves  and  in  their  simplest  form  are  referred  to  as  a  Mach 
stem. 

The  other  major  result  is  the  small  amount  of  computer  time  required  for  these  simulations.  The  waves 

are  significantly  nonlinear  with  the  ratio  of  significant  wave  height  to  water  depth  given  by 

Hs/h  =  l.5m/8m  =  0.lS75  mdthe  maximum  wave  height  to  depth  is  # max  l  h  =  3.4m  /  8  m  =  0.425 

The  integration  domain  is  500  m  by  500  m  and  has  128  x  128  spatial  bins.  A  total  of  500  time  values 
were  computed  to  make  the  film  (Til  send  you  a  .avi  file  if  you  request  it).  The  total  cpu  time  was  1 
min  10  sec.  This  compares  to  a  split-step  FFT  code  which  took  4.5  hours  to  compute  the  same  problem. 
Thus  the  preprocessor  multi-dimensional  Fourier  model  which  I  have  developed  is  about  230  times 
faster  than  the  split-step  FFT  code.  For  a  number  of  reasons  which  I  will  not  mention  here,  I  anticipate 
another  order  of  magnitude  improvement  in  computer  time  in  the  preprocessor  part  of  the  model  over 
the  next  year  or  so. 


Fig.  2  Initial  conditions  for  simulation  of  KP  equation.  The  significant  wave  height  is 

1.5  m  and  the  water  depth  is  8  m. 
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Fig.  3  A  large  rogue  wave  of  height  3.4  m  in  the  Boussinesq  (xKP)  simulation.  The 
significant  wave  height  is  1.5  m.  The  rogue  waves  have  maxima  shown  in  red. 


Fig.  4  Another  large  rogue  wave  of  height  3.3  m  in  the  Boussinesq  (xKP)  simulation. 

The  significant  wave  height  is  1.5  m. 
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RESULTS 


The  numerical  model  of  the  Euler  equations  uses  particular  transformations  similar  to  (2)  plus 
additional  other  mathematical  machinery  to  simulate  the  wave  motion  simultaneously  in  terms  of  both 
the  surface  elevation  and  the  velocity  potential.  There  are  three  theta  functions  that  have  to  be 
accounted  for  and  collapsed  onto  the  linear  Fourier  modes.  One  of  the  most  important  considerations  is 
that  the  theta  function  now  has  a  three  dimensional  form: 

...  .  _  \  '  /m-kv.v-/m-ki  r+/m-U  r-/ni-(i)  /  ■  /m-<j) 

me 

me0 

where  the  “nomes”  are  given  by: 

'  m- Bin 

=  e  2 

But  in  the  formulation  we  need  to  evaluate  this  function  at  the  free  surface,  z  -  _  This  gives  for 

the  theta  function  at  the  free  surface : 

0(x,y,T),t)  =  I  ^  Jm'kxx+im'\Lyy+im»\Lz7j(x,y,t)-im'(o  /+zm*(|) 
me0 

It  is  tempting  to  write: 

eim.kzjj(x,y,t)  =  j  +  -  I(m.kz)2  t]2(x,y,t)  -  -^(m.kz)3  t]3(x,y,t )  + ... 

2  6 

which  should  converge  provided  that  the  wave  steepness  has  the  property 

m.kz?](x,y,t)«l 

This  can  of  course  be  substituted  into  the  theta  function  to  obtain 
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0(x,y,z,t )  =  ■Z  ^  ^/m-kvx+/m-kv3;+/m-kz77(x,3;,/t)-/m-cc)  £+/m-<|>  ^ 

me0 

m<E0 

me0 


-^3(x,,,0Z(m-k,)3^m- 

me0 


kv.v-/m-k1.r+/m-k./7(.v,r,/  )-/m-co  /  ■  /m-(j> 


+  ... 


These  are  the  fundamental  types  of  multi-dimensional  Fourier  series  which  appear  in  the  formulation. 

A  flowchart  of  the  Euler  equation  simulation  is  shown  in  Fig.  5.  In  Fig.  6  we  show  the  appearance  of  a 
rogue  wave  from  one  of  the  simulations  in  deep  water. 

IMPACT/APPLICATION 


The  impact  of  this  research  will  occur  in  general  for  the  nonlinear  Fourier  analysis  of  shallow  and 
deep-water  wave  trains.  Specific  results  will  provide  for  a  deeper  understanding  of  nonlinear  wave 
dynamics. 
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Solving  Nonlinear  PDEs  With 
The  Multi-dimensional  Fourier 
Transform 


Fig.  5  Simple  flow  chart  of  the  numerical  model  for  the  Euler  Equations. 


Fig.  6  A  single-peak  rogue  wave  with  deep  holes  on  the  leading  and  trailing  edges.  The 
vertical  scale  should  be  multiplied  by  10.  The  height  of  the  wave  is  46  m  and  has  a  wave 

length  of 320  m. 


TRANSITIONS 

Transitions  expected  are  related  to  the  use  of  the  codes  as  guidance  to  ships  and  unmanned,  unteathered 
vehicles  as  the  kind  of  environment  in  which  one  resides  and  for  the  real  time  sampling  of  the 
environment,  including  the  acoustic  environment. 

RELATED  PROJECTS 

An  intimate  relationship  between  our  results  and  other  projects  exists  because  the  sea  surface  provides 
a  major  forcing  input  to  many  kinds  of  offshore  activities,  including  the  dynamics  of  floating  and 
drilling  vessels,  barges,  risers  and  tethered  vehicles.  The  present  work  leads  to  a  nonlinear 
representation  of  the  sea  surface  forcing  and  vessel  response  for  shallow  water  waves. 
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